Vibronic effects in single molecule conductance: First-principles 
description and application to benezenealkanethiolates between 

gold electrodes 

C. Benesch, 1 M. Cizek, 2 J. Klimes, 2 M. Thoss, 1 and W. Domcke 1 

1 Department of Chemistry, Technical University of Munich, 

Lichtenbergstr. 4, D-85747 Garching, Germany 
2 Charles University, Faculty of Mathematics and Physics, 
Institute of Theoretical Physics, Prague, Czech Republic 

Abstract 

The effect of vibrational motion on resonant charge transport through single molecule junctions 
is investigated. The study is based on a combination of first-principles electronic structure calcu- 
lations to characterize the system and inelastic scattering theory to calculate transport properties. 
The extension of the methodology to describe hole transport through occupied molecular orbitals 
is discussed. The methodology is applied to molecular junctions where a benzene molecule is con- 
nected via alkanethiolate bridges to two gold electrodes. The results demonstrate that, depending 
on the coupling between the electronic 7r-system of the benzene ring and the gold electrodes, 
vibronic coupling may have a significant influence on the transport properties of the molecular 
junction. 
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I. INTRODUCTION 



Recent advances in experimental studies of single molecule conduction [l, 2, ^, Q, O, E , 
9! have stimulated great interest in the basic mechanisms which govern electron transport 
through nanoscale molecular junctions [lO, 
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12|. An interesting aspect that distinguishes 



molecular conductors from mesoscopic semiconductor devices is the possible influence of the 
nuclear degrees of freedom of the molecular bridge on electron transport. Due to the small 
size of molecules, the charging of the molecular bridge is often accompanied by significant 
changes of the nuclear geometry. The current-induced excitation of the vibrations of the 
molecule may result in heating of the molecular bridge and possibly breakage of the junction. 
Conformational changes of the geometry of the conducting molecule are possible mechanisms 
for switching behavior and negative differential resistance fl3| . Furthermore, the observation 
of vibrational structures in conduction measurements allows the unambiguous identification 
of the molecular character of the current. 

Vibrational structures in molecular conductance were observed, for example, in electron 

n 

transport experiments on H 2 , HD, and D 2 between platinum electrodes [14| using the me- 
chanically controllable break junction technique. Thereby, changes in the conductance were 
assigned to switching between two different local geometrical configurations induced by the 
transverse motion of the molecule. The observed vibrational structures could be classified 
as longitudinal or transversal 15]. In studies of C 60 molecules between gold electrodes 
the center of mass motion of the molecular bridge was observed. Other experiments on 



this system If], [l7| as well as on C 70 |7|, (070)2 [l8j], and copper phtalocyanin [6| on an 
aluminum oxide film showed structures which were related to the internal vibrational modes 
of the molecule. The aluminum oxide layer in the latter experiment acted as an insulating 
layer, which effectively reduces the electronic coupling between the molecule and the metal 
substrate thus facilitating the effect of vibrational motion on the conductance. In another 
STM experiment 19] , pronounced progressions of vibrational modes were observed and the 
dependence of the vibronic coupling on the spatial position of the STM-tip was demon- 
strated. Moreover, vibrational signatures of molecular bridges have also been observed in 
off-resonant inelastic electron tunneling spectroscopy (IETS) {20! . 21]. 

This tremendous experimental progress has inspired great interest in the theoretical 
modelling and simulation of vibrationally-coupled electron transport in molecular junc- 
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371 ] . In the low- voltage, 



off-resonant transport regime combinations of electronic-structure calculations and non 
equilibrium Green's function theory (NEGF), employing the self-consistent Born approx 
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imation (SCBA), have been used to investigate vibrational signatures in IETS 
Other approaches used to describe inelastic features in the off-resonant region include scat- 
tering theory combined with an expansion of either the molecular Green's function to first 



order in vibronic interaction 



28j | or the systems wavefunction in the nuclear coordinates 



29]. 



The majority of the studies of vibronic effects in the resonant tunneling regime (for higher 
voltages) have been based on generic tight-binding models, using an NEGF approach based 
on the equation of motion method on the Keldysh contour 30j , or kinetic rate equations to 



calculate the current 
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361 ] . These model studies have demonstrated that 



the vibrational motion of the molecular bridge may affect the current-voltage characteristics 
significantly. 

Many of the studies of vibrationally-coupled electron transport reported so far invoke 
approximations which restrict their applicability either to small electronic-vibrational cou- 
pling, small molecule-lead coupling or separability of the vibrational modes. To circumvent 
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401 ]. which is based on inelastic 



this limitation, we recently proposed an approach 
multi-channel scattering theory and the projection-operator formalism of resonant electron- 
molecule scattering [41]. Within the single-electron description of electron conduction, this 
approach is valid for strong electronic-vibrational and molecule-lead coupling and thus allows 
the study of electron transport in the resonant regime. 



In previous work 38|, ]39j], we have applied this formalism to study vibronic effects in 



molecular conductance based on generic models for molecular junctions. In a recent short 
communication [42], we have combined the scattering theory approach with first-principles 
electronic-structure calculations to characterize the molecular junction. In this paper, we 
give a detailed description of the methodology including a derivation of the model Hamil- 
tonian for hole transport, a description of the partitioning method used to determine the 
parameters of the Hamiltonian by first-principles electronic structure calculations and an 
outline of the scattering theory approach. The methodology is applied to study vibronic 
effects on conductance in molecular junctions where a benzene ring is connected via alka- 
nethiol bridges to two gold electrodes (Au-S-(CH 2 ) n -C 6 H 4 -(CH 2 ) n -S-Au, cf. Fig. HJ. This 
class of molecular junctions was chosen, because the change of the length of the alkyl chain, 



3 



n, allows a systematic variation of the coupling between the electronic 7r-system of the 
phenyl ring and the gold electrodes and thus a variation of the lifetime of the electron on 
the molecular bridge. Our previous studies show that, due to the strong molecule lead 

ate (corresponding to n — 0) between gold electrodes exhibits only 



coupling, benzenethio 

minor vibronic effects 42j. In this work, we consider benzene-di(ethanethiolate) (n = 2) and 
benzene-di(butanethiolate) (n = 4), which correspond to moderate and weak molecule-lead 
coupling, respectively. 

II. THEORY 

A. Charge transport model Hamiltonian 

To study charge transport through single molecule junctions employing inelastic scat- 
tering theory, we use an effective single-particle model Hamiltonian with parameters deter- 
mined by first-principles electronic structure calculations. In this section we briefly outline 
the derivation of the Hamiltonian. 

The Hamiltonian of a metal- molecule- metal (MMM) junction including the electronic- 
nuclear interaction is given by the generic expression (we use atomic units where h = e = 1), 

H = T n + V nn (R)+H e (R). (2.1) 

Here T n denotes the kinetic energy of the nuclei, 

T n = V" -EL (2.2) 

a 

where p a is the momentum and M a is the mass of nucleus a. The Coulomb interaction of 
the nuclei is given by 

where Z a is the atomic number and R a denotes the Cartesian coordinates of the position of 
nucleus a. The electronic Hamiltonian, if e (R) includes the kinetic energy and the Coulomb 
interaction of the electrons as well as the electronic-nuclear Coulomb interaction and depends 
parametrically on the nuclear coordinates R. 

We first focus on the electronic degrees of freedom and consider the equilibrium situation 
without bias voltage at zero temperature. In this case, the junction is in its electronic ground 



state. Employing an effective one particle description, the ground state is given by a single 
Slater determinant, \^f g ) = n!Si c l|0}> where |0) denotes the vacuum state, the operator 
c] creates an electron in the single particle state (molecular orbital) i and the product is 
taken over all n e i electrons in the neutral junction. \ty g ) is an approximate solution of the 
electronic Schrodinger equation, 

H e (R)\* g ) = E g (R)\S> g ) (2.4) 

with eigenenergy E g (R), which parametrically depends on the nuclear coordinates. The 
energy E g (R) + V nn (R) represents the adiabatic potential energy surface of the electronic 
ground state, which will serve as a reference state in the following. 

If an external voltage is applied to the MMM-junction, electrons (or holes) are transferred 
from an external reservoir to the junction and vice versa, i.e. the number of electrons on 
the junction, n e i, will change. Here we consider charge transfer processes which result in 
one additional electron (corresponding to the anion) or one additional hole (corresponding 
to the cation) on the junction. Within Koopmans' theorem the energy of a cation is given 
by subtracting the energy of the single-particle state (molecular orbital), from which the 
electron was removed, from the energy of the neutral ground state. The energy of an anion, 
on the other hand, is given by adding the energy of the single-particle state that was occupied 
to the energy of the neutral ground state. An electronic Hamiltonian that can describe both 



situations is given by [43] 



H e (R) = E g (R)- E AR)cjc}+ E i( R ) c k, (2-5) 

j'Gocc. iSunocc. 

where the index j labels all orbitals that are occupied in the neutral reference state, while 
the index % labels all orbitals that are unoccupied in the neutral reference state. Thereby, 



effects due to electron-electron correlation and orbital relaxation have been neglected [441 ]. 

If the electronic spectrum of a molecular junction is such that the levels, which contribute 
to conduction, are all located far below or far above the Fermi energy, the third or second 
term in Eq. (I2.5P can be neglected, respectively. In the latter case the current is dominated 
by electron transport, while in the former case hole transport prevails. Thus in the two 
limiting cases, we obtain the electronic Hamiltonian 

H e (R) = E 9 (R) + Ei{R)c\ Ci , (2.6a) 



H e (R) = E g (R) - EAB)c jC ). (2.6b) 

3 

for electron or hole transport, respectively. In all applications considered below, the trans- 
port is dominated by hole transport. Therefore, we will consider in the following only this 
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case. A more detailed discussion and comparison of both cases is given in Ref. 

To describe charge transport trough a molecular junction it is expedient to partition the 
overall system into the molecule and the left and right leads. As will be described in more 
detail in Sec. IIIB[ this results in a splitting of the second term in Eq. (12.6bl) into three terms, 
describing the molecule, the left electrode, and the right electrode separately, 



fle(R) = E g (R)- E m {R)c m cl- E kCkcl 

m£M k£L,R 
Y {yrnkC k C ] m + V km C m c[^ . (2.7) 



E 

mSM k£L,R 



Here the indices m and k denote orbitals that are either localized on the molecule (M), or 
on the left (L) or right (R) lead, respectively, and and c\ are the corresponding creation 
operators, and V mk describes the coupling between molecule and leads. In Eq. ( 12. 7L we 
have assumed that direct interactions between left and right lead can be neglected and 
that the energies of the lead states E^, as well as the coupling matrix elements, V m k, are 



46| 



approximately independent on the nuclear coordinates 

In the applications considered below, we use the harmonic approximation for the nuclei. 
In this approximation, the potential energy of the nuclei in the electronic ground state, 
E g (R) + V nn (R), is expanded up to second order in the nuclear coordinates around the 
equilibrium geometry (R eg ) of the neutral molecule. This results in the following expression 
for the nuclear Hamiltonian of the neutral reference state 

T n + V nn {R) + E g {R) = H n0 = J2 [<4 a i + I J > ( 2 - 8 ) 

where u>i is the frequency of the Zth vibrational normal mode of the molecule with dimension- 
less coordinate qi and a',, ai denote the corresponding creation and annihilation operators. 
Employing, furthermore, an expansion of the molecular energies Ej(R) up to first order 
around the equilibrium geometry of the neutral molecule, the overall Hamiltonian for hole 
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transport is given by 

Am) 



( 1\ K (m> ( \ 

H = Yl Ul ( °i flj + 2 ) ~ E ™ C ™ C m - ^ -W + a «) C ™4 

J ^ ' meM l,m€M * 

- ^ E kCkc{ + ^ {yrnkC ] m C k + Vjfc m 4c m J . (2.9) 

keL,R meM keL,R 

Here, we have introduced the electronic-vibrational (vibronic) coupling parameters 

The linear form of the vibronic coupling implies that only normal modes belonging to the 
totally symmetric representation of the respective symmetry group have a non-vanishing 
coupling constant. The harmonic approximation of the nuclear potential employed in Eq. 
(12. 9p is valid for small amplitude motion around the equilibrium geometry. For a treatment 
of large amplitude (e.g. torsional) motion in molecular junctions, see Ref. [39. 



The Hamiltonian derived above can be used within different methods to study transport 

n 

properties of molecular junctions, e.g. density matrix theory [47J or nonequilibrium Green's 
function approaches 48] . In this paper we use it in combination with inelastic scattering 
theory. When employing scattering theory it is expedient to represent the Hamiltonian in 
terms of effective single-particle states for the electronic degrees of freedom. In the case 
of hole transport, the relevant single particle states are \(/> m ) = c m \^/ g ), \4>k) = Ck\^ g ), 
which describe an additional hole on the molecule or the leads, respectively. The effective 
single-particle representation of the Hamiltonian is obtained by projecting Eq. (12.91) onto 
the single-particle space, i.e. H — > Yin' 1 0i) (0*1-^1 I This results in 



H = H n0 — H M — H ne — H L — H R — V (2-11) 



with 
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5>i(ata + ^)> (2.12a) 



H M = J2 \<Pm)E m (<p m \, (2.12b) 

meM 

(m) 

H ™ = ^/f{ ai + V \M(^ m \, (2.12c) 



LmeM 



H L + H R = ^ \<j> k )E k (<f> k \, (2.12d) 



keL,R 



V = J2 (l^Vkmi^kl + \(pk)V mk ((p m \) ■ (2.12e) 



meM keL,R 



The negative signs in Eq. (12. lip are due to the fact that we consider hole transport. The 
Hamiltonian for electron transport can be derived in a similar way and is explicitly given in 
It is noted, that for purely electronic transport calculations (i.e. without vibronic 



Ref. 
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coupling, nf 11 ^ = 0), the different Hamiltonians for electron and hole transport give the same 
current-voltage characteristic. Including vibronic coupling, this is, however, no longer the 
case 45] . 



B. First-principles determination of parameters 

The parameters of the Hamiltonian (12. lip describing the molecular junction are deter- 
mined by first-principles electronic structure calculations. To this end, the overall junction 
is partitioned into five parts: The molecular part, the part of the left and right electrodes 
that is explicitly included in the quantum chemistry calculation (in the following referred to 
as the contacts), and the remaining parts of the left and right electrodes. The molecule and 
the contacts form the extended molecule and are treated explicitly by quantum chemistry 
calculations. The influence of the remaining part of the electrodes is taken into account 
implicitly employing the surface self energy of a gold (111) surface. 



1. Electronic structure methods 



The parameters for the explicitly calculated part of the molecular junctions were deter- 
mined employing electronic-structure calculations performed with the TURBOMOLE pack- 
age (V5-7) [49f using density functional theory (DFT) with the B3-LYP hybrid functional 
and the SV(P) basis set including ECP-60-MWB on the gold atoms. 

The final structure of the three systems considered (cf. Fig. [1]) was the result of several 
optimization cycles: A full geometry optimization of the isolated neutral molecule was fol- 
lowed by a full geometry optimization of the system obtained after replacing the hydrogens 
on sulfur by two gold atoms. Covalent bonding to t wo g old atoms is the preferred bond 
formation if no symmetry constraints are applied {4], 50]. Finally, a second metal layer 
containing five gold atoms was added on both sides and the geometry of the system was 
again optimized, thereby keeping the internal coordinates of the second gold layer fixed. 
The final geometry of the metal contacts was obtained by increasing the second layer by 



S 



several atoms and/or adding two more metal layers, which are cutouts of the (lll)-plain of 
the face-centered-cubic lattice of solid gold. A single point DFT calculation was performed 
on this system. The Fermi energy was approximated by the average value of the energies of 
the HOMO and LUMO. 



2. Partitioning method and electronic parameters 

The electronic structure calculations discussed above result in delocalized molecular or- 
bitals for the extended molecule. To determine the single-particle states \<fij), \4>k) employed 
in the Hamiltonian (12. lip , which are localized on the leads or the molecule, respectively, 
and the corresponding electronic energies and coupling parameters, we employ a partition- 
ing method. This method is based on the following four steps: (i) separation of the overall 
Hilbert space into left contact, molecule, and right contact, (ii) partitioning of the Hamilto- 
nian according to this separation, (iii) introduction of a self energy to describe the influence 
of the remaining infinite electrodes (iv) separate diagonalization of the three blocks of the 
partitioned Hamiltonian. The detailed procedure is discussed in the following. 

As was mentioned above, in the present paper we work within an effective single-electron 
picture. Thus we identify the effective electronic Hamiltonian with the Kohn-Sham matrix, 
F, of the self-consistent-field DFT calculation. The Kohn-Sham matrix is represented in 
a set of atomic orbitals \Xn)- To separate the overall Hilbert space into three mutually 
orthogonal subspaces describing the molecule and the two contacts, the following projection 
procedure was employed: First, the set of atomic orbitals of the overall system, {|Xn)} is 



divided into three groups, which are centered at the left contact, the molecule, and t 



contact, respectively. Since it is advantageous to work with orthogonal orbitals 
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re right 



set of atomic orbitals of the overall system is orthogonalized according to Lowdin 
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\x 



i> = J2( S ~ 1/2 )"i\Xi) , (2-13) 



where S denotes the atomic orbital overlap matrix with elements Si n = (xi\Xn)- The new 
basis functions obtained, \xn), exhibit a minimal deviation from the original ones in a least- 
square sense and hence, their localization is preserved. In particular, the classification into 
the three groups (left contact, molecule, right contact) is still valid. 

The new set of orthogonal basis functions {|Xn)} is then used to partition the Kohn-Sham 
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matrix from the converged DFT calculation into the three subspaces. The Kohn-Sham 
matrix in the orthogonal basis is given by 



F = S-2FS-*. 



(2.14) 



According to the separation of the orthogonalized basis set to the three groups - left con- 
tact, molecule, right contact - the transformed Kohn-Sham matrix can be arranged in the 
following block structure 



^ Fl Flm Flr ^ 



V 



(2.15) 



/ 



Fml Fm Fmr 
Frl Frm Fr 

The thus obtained Kohn-Sham matrix describes the extended molecule, which includes a 
part of the electrodes. To incorporate the influence of the remaining (infinite) parts of the 
electrodes, we use a surface Green's function technique. Within this method, the influence of 
the remaining part of the electrodes is described by the corresponding self energy. Within the 
Green's-function formalism employed below to calculate transport properties this amounts 
to adding the self energy to the Kohn-Sham matrix elements belonging to the gold atoms 
of the outermost layer of the contacts 



rrtot 



F + S sf . 



(2.16) 



The surface self energy S sf was obtained by calculating the surface Green's function for a 
(111) gold surface representing the semi-infinite electrode |55j. The surface Green's function 
of a gold (HI) surface was determined using a generalization of the method developed by 



561 ] . Thereby, the (111) surface of gold was described by a tight binding model 



Sancho et al. 

[57I . Thereby, we have neglected nondiagonal contributions, which correspond to coupling 
between different gold atoms. It is noted that due to the added self energy, the Kohn-Sham 
matrix is no longer real but complex symmetric 

In the final step, the three different blocks of the Kohn-Sham matrix F tot are diagonalized 
separately. Denoting the matrices of right eigenvectors for the three subspaces as Ul, Um, 
Ur, respectively, the Kohn-Sham matrix is transformed to the new basis of eigenvectors 



F 



U T F tot U 



^ El V L m Vlr ^ 

Vml Em Vmr 
Vrl Vrm E r 



(2.17) 
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where 



' U L ^ 



V 



/ 



V = U M • (2.18) 
o U R 

The elements of the thus transformed Kohn-Sham matrix F determine the electronic pa- 
rameters of the model Hamiltonian (12.111) : The diagonal submatrices El and Er contain 
the lead energies E^, the diagonal submatrix Em contains the eigenenergies of the molecular 
bridge and the submatrices V describe the electronic coupling matrix elements Vk m between 
states localized on the leads and the molecule. In all our model systems, the distance be- 
tween the two metal clusters is relatively large, Therefore, direct lead-lead coupling was 
small enough to be neglected. 

The electronic states in the Hamiltonian (12.111) can be identified with the eigenstates of 
the three different blocks of the Kohn-Sham matrix and are related to the original atomic 
orbital basis via 



*) = J2 u ™( s ~ l/2 Uxi), 



(2.19) 



where a G L, M, R. While the eigenvectors of F\ 



tot 

M > IVm 



) , form an orthogonal basis in molec- 



ular space, the eigenvectors of the non-Hermitian matrix F^, \<j>k), form a bi-orthogonal 



basis in the left lead space. Their overlap, 



/) = S^Skk', where \4>k) is the right eigen- 



vector of _F£ ot and 



. | is its transposed, is in general a complex number. The same is true 
for the right lead space. We note that a similar separation scheme was employed in the 
study of electron transfer dynamics in dye-semiconductor systems 581 ]. 



3. Nuclear parameters 

To characterize the nuclear degrees of freedom of the molecular bridge and determine 
the corresponding parameters in the Hamiltonian (12.111) a normal mode analysis on an 
extended molecule including two gold layers (i.e. five gold atoms on each side) was performed. 
Thereby, the atomic mass of the gold atoms was set to 10 9 atomic mass units to separate 
the nuclear motion in the molecular bridge from that of the gold clusters. 

The electronic-nuclear coupling constants were obtained from the numerical gradi- 
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ents of the energies Ej with respect to the dimensionless normal coordinates qi, 

M = E m (+A qi ) - E m (-A gi ) 

1 2A Ql ' 1 ' ' 

To this end, two DFT calculations with molecular geometries elongated by ±A<# = ±0.1 
from the equilibrium geometry were performed. The value of Aqi = 0.1 was chosen based 
on convergence tests. 



C. Observables of Interest 



Several observables are of interest to characterize the effects of vibrational motion on 
charge transport through single molecule junctions. Here, we will focus on the inelastic 
transmission probability, which describes the transmission process of a single electron or 
hole through the molecular junction, and the current- voltage characteristics. In all the 
examples considered below, the transport is dominated by hole transport. The inelastic 
transmission probability of a single hole from the left to the right lead as a function of initial 
and final energy of the hole is given by the expression 

tw^, E f ) = 4tt 2 J2 E E p ^ s ( E »f ~ E f- E »< + E i) 

x5(E t - E h )S(E f - E k} ) \{u f \{k } \VG{E)V\h)\u^ . (2.21) 

Here, the 5-function accounts for energy conservation, E = E Vi — Ei = E Uf — Ef with E Ui , 
and E Uf being the energy of the initial and final vibrational states respectively, 
—Ei and — Ef are the initial and final energies of the hole, G(E) = (E + — if) -1 is the 
Green's function, and P Vi = (vi\po\vi) denotes the population probability of the initial 
vibrational state, p = e ~ Hn0 ^ kBT y Z . In the systems considered in this work, the hole 
couples primarily to modes with relatively high frequencies. As a result, thermal effects 
are not expected to be of relevance and the initial vibrational state is assumed to be the 
ground state, i.e. po — |0)(0|. The total transmission probability, Tn^L^Ej), is obtained by 
integrating TR < _i(£ , i , Ef) over the total range of final energies of the hole. 

T R ^ L {Ei) = J dE f T R ^ L (E t ,E f ) . (2.22) 

To calculate the transmission probability, the Green's function G(E) in the expression 
for the transmission probability is projected onto the molecular space using the projection 

12 



operators P = E mG Af 1 0m) (0m I, Ql = EfceL \<t>k) {<f>k\, Qr = T,keR I0fc)<0fcl> as wel1 as the 
Lippmann-Schwinger equation G = Go + GoV^G. This results in the expression 

T R ^ L (E U E f ) = p vAEu t -E f - E Vi + E t ) (2.23) 
xtr M {(vi\T L (-Ei) G\ fl {E v% - Ei) \u f ){u f \ T R (-E f ) G M (E Vi - ^)l^>} , 



where the trace is taken over the electronic states on the molecule (M) and the Green's 
function projected onto the molecular bridge is given by 

G M (E) = PG(E)P (2.24) 
= [E+ - PHP - - H n0 ) - ££ ad (£ - H n0 )] 1 

= [E + + H M — H n0 + H ne — T^(E — H n0 ) — T}ft d (E — H n0 )] 

Here, ^[f R (E) denote the self energy due to coupling to the left and right lead, respectively 

£^ ad (£) = PVQ L (E + + H L y l Q L VP = - l -T L {E) + A L (E), (2.25a) 
££ ad (£) = PVQ R (E + + H R )- l Q R VP = - l -T R {E) + A R (E), (2.25b) 

where Tl(E) and Al(E) are the corresponding width and level-shift functions. The matrix 
elements of the self energy in the molecular space are given by 

k 

It should be emphasized that the thus defined self energy S 1 £ ad (i?) corresponds to hole 
transport. It is related to the commonly used self energy for electron transport, H l ^(E) 
via 

(^r d (E))mn = - (EKMOW ■ (2-27) 

In the practical calculations, the latter self energy is obtained using the partitioning method 
outlined above and is explicitely given by 

(PtfWU = £ S U^-E k ) (2 ' 28) 

and analogously for S^i?. The appearance of the overlap factor S% in Eq. (I2.28P is a 
consequence of the bi-orthogonal basis on the leads (cf. the discussion in Sec. Ill B 2j) . 
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It is noted, that in contrast to purely electronic transport calculations as well as appli- 



cations of the non-equilibrium Green's function formalism to vibronic transport [25|, [59], 
the Green's function Gm and the self energies S^ ad , £^ ad in Eq. (12.25b[) are operators with 
respect to both the electronic and nuclear degrees of freedom. Thus, the Green's function 
has to be evaluated in the combined electronic-vibrational Hilbert space. In the calculations 
reported below, the Green's function Gm is obtained by inverting a basis representation of 
the operator E + — PHP — S^ ad (i? — H n0 ) — S^ ad (i? — H n0 ) for each energy E employing a 
harmonic oscillator basis for the vibrational modes. 

Based on the transmission probability (Eq. 12.231) . the current through the molecular junc- 



tion is obtained using a generalized Landauer formula [60] 



2e 
~h 



j dE t j dE f {T R ^ L (E t ,E f )f R (E f )[l - f L (E t )] 

- T L ^ R {E h E f )f L (E f )[l - f R (Ei)]} , (2.29) 



where fi(E), f R (E) denote the Fermi distribution for the left and the right lead, respectively, 
taken at zero temperature. 

While equations (12.211) and (12.231) for the single-hole transmission probability involve 
no approximation, expression f!2.29j) for the current is valid if many-electron processes are 
negligible. In particular, non-equilibrium effects in the leads and electron correlation due 
to electronic-vibrational coupling are not taken into account. Furthermore, it is implicitly 
assumed, that the nuclear degrees of freedom of the molecular bridge relax to the vibra- 
tional equilibrium state after scattering of a hole. For the applications considered below, 
comparison with calculations based on nonequilibrium Green's function methods indicate 



that these assumptions are justified 48|. Without vibronic coupling the expression for the 



current, Eq. (I2.29p . reduces to the usual Landauer formula [61]. 

In principle, the basis states \<fi m ), \4>k), the hole energies and the nuclear parameters 
depend on the bias voltage. For the studies in this work, we did not include the voltage 
self-consistently in the DFT calculation but used, for simplicity, parameters obtained from 
electronic-structure calculations at equilibrium and assumed that the bias voltage V enters 
the formulas only via the chemical potentials of the leads Hl/r = e / =t e ^/2, where €f 
denotes the Fermi energy. The energies of the lead states for finite voltage are thus given by 
Ek ± eV/2. Since we do not invoke the wide-band approximation, the Green's function Gm, 
the self energies ^ l [J R as well as the width functions T L / R also depend on the bias voltage. 
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FIG. 1: The three molecular junctions investigated: benzene-di(ethanethiolate) (BDET) between 
pyramidal gold contacts (top), BDET between cuboid-shaped gold contacts (middle), and benzene- 
di(butanethiolate) (BDBT) between pyramidal gold contacts (bottom). 

III. RESULTS AND DISCUSSION 

The theoretical methodology outlined above has been applied to study vibrational ef- 
fects in electron transport through the three different molecular junctions depicted in Fig. 
[U which include benzene-di(ethanethiolate) (BDET) between pyramidal and cuboid-shaped 
gold contacts as well as benzene-di(butanethiolate) (BDBT) with a pyramidal contact ge- 
ometry. To reduce the computational effort, in the calculations presented below, the four 
vibrational modes with the strongest vibronic coupling (as determined by the ratio of the 
electronic-vibrational coupling k and the electronic coupling T) were explicitly taken into 
account. Furthermore, the number of electronic states \<p m ) on the molecular bridge, which 
were explicitly included in the calculation, was reduced by including only those with energies 
in the vicinity of the Fermi energy (12 for BDET and 6 for BDBT). 
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FIG. 2: Normal modes of BDET included in the calculation. 

TABLE 1: Parameters of the four most important vibrational modes of BDET between pyramidal 
gold contacts including frequencies, periods, and vibronic coupling in the electronic states A and 
B. 





uj (cm 1 ) 


T(fs) 




K (B) ( meV ) 


(a) 


544.48 


61 


76 


22 


(b) 


1197.11 


28 


51 


69 


(c) 


1229.49 


27 


110 


47 


(d) 


1671.56 


20 


136 


162 



A. Benzenedi(ethanethiolate) with a pyramidal gold cluster geometry 

We first consider the transport characteristics of BDET between two pyramidal gold clus- 
ters (Fig.fU top panel). Fig. [2] depicts the four vibrational normal modes of BDET with 
the strongest vibronic coupling included in the calculation, which can be characterized as: 
C-C-C bending ((a) and (c)), C-C-H bending (b), and C-C stretching (d). The correspond- 
ing normal mode frequencies and vibronic coupling constants in the two most important 
molecular orbitals are given in Tab. [TJ 

The transmission probability of BDET at zero voltage is depicted in Fig. [3] (a). In addition 
to the transmission probability based on a vibronic calculation, also the result of a purely 
electronic calculation (where all vibronic coupling constants k\ have been set to zero) is 
shown. The result of the purely electronic calculation exhibits two pronounced resonance 
peaks at energies -2.23 eV and -1.85 eV. These two resonances are the closest to the Fermi 
energy and thus determine the transport process. The two resonances are caused by the 
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FIG. 3: (a) Total transmission probability through a BDET molecular junction with a pyramidal 
gold cluster geometry at zero voltage as a function of the initial energy of the hole (relative to the 
Fermi energy). The two orbitals, denoted A and B, dominate the transmittance at the indicated 
peaks. Only the energy range with non-zero transmittance in the interval [-2.5:+2.5], corresponding 
to a voltage window of 5 V, is shown, (b) Current-voltage characteristic of BDET bound to 
pyramidal gold contacts. Shown are results of calculations with (solid line) and without (dashed 
line) molecular vibrations. Only the non-zero part of the positive voltage range is shown. The 
different lines in the two panels show results of calculations with (dashed line) and without (solid 
line) vibronic coupling. 

orbitals A and B depicted in Fig. [3] (a). Assuming that the projected molecular space 
includes the same number of electron pairs as in the isolated molecule, the states A and 
B would correspond to the HOMO-4 and HOMO-5 of the isolated molecule. While orbital 
A resembles an ei s -orbital of benzene, orbital B has additional contributions on the ethyl- 
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groups and the sulfur atoms. As a consequence, state A has smaller electronic coupling 
to the leads (Taa = 1.9 • 10~ 4 eV, corresponding to a lifetime of r = Ti/Taa = 3466 fs) 
resulting in a very narrow peak in the transmission probability, whereas the significant 
coupling of orbital B to the leads (Tbb — 1-1 ■ 10 _1 eV, corresponding to a lifetime of r = 6 
fs) results in a rather broad structure. It is noted that in this estimate of the lifetime, 
the width function T has been taken at the corresponding resonance peak energy. Due to 
the interaction with the leads and neighboring states, the peak energy, in general, does 
not coincide with eigenenergies of the states after projection, E m . For example, there exist 
four pairwise degenerate sulfur p-states with energies close to the Fermi level, which do not 
appear themselves in the transmission probability but influence the transmission via their 
interaction with states A and B. This interaction is mediated by the gold clusters via the 
non-diagonal elements in the self energy matrix. 

The comparison between the results of vibronic (solid line in Fig. [3] (a)) and purely elec- 
tronic (dashed line) calculations demonstrates that the electronic-vibrational coupling in 
BDET alters the transmission probability noticeably. In particular, it results in a shift of the 
electronic resonance peaks due to nuclear relaxation and the appearance of new structures, 
which correspond to vibrational states in the molecular cation. Specifically, the electronic 
resonance of state B is shifted by 0.08 eV to the peak corresponding to the 0° transition 62] 
and the transmission shows resonance structures at -1.92 eV and -1.98 eV corresponding to 
single vibrational excitations of mode (b) {{vh)\ transition) and (d) ((^ c )o transition). It is 
noted that in contrast to previous calculations 42j, where charge transport through BDET 
was described as electron instead of hole transport, the shift due to nuclear relaxation is 
in the opposite direction and the higher excited vibrational states in the vibrational pro- 
gressions are located at lower energies. The 0q transition peak belonging to state A, which 
may be expected to be a narrow line, cannot be seen because it is coincidentally shifted 
into resonance with a sulfur p-state transmission peak. This peak, which does not appear 
directly in the transmission probability, causes a broadening of the resonance of state A. 
The structure resulting from this effect can be seen around peak A but it is too broad to 
assign individual modes. 

The current-voltage characteristic of BDET is shown in Fig.|3] (b). Due to the symmetry 
of the junctions, the current fulfills the relation J(— V) = — J(V) and, therefore, only the 
positive voltage range is shown. The result obtained from a purely electronic calculation 
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FIG. 4: Conductance of BDET based on a purely electronic (dashed line) and a vibronic (solid 
line) calculation. 



(dashed line) shows an increase of the current at about 3.5 V related to hole transport 
through state B. The increase is followed by a pronounced decrease of the current at 4 V. 
This negative-differential resistance (NDR) effect is a result of the voltage dependence of the 
self energies S^ ad , S^ ad and the corresponding width functions T and will be discussed in 
more detail below. The weakly coupled state A results only in a small step-like increase of 
the current shortly before 4.5 V and the following rise of the current is due to the increasing 
width of the resonance peak due to state B which appears at -1.85 eV in the transmission 
probability. 

Including the coupling to the nuclear degrees of freedom changes the current-voltage 
characteristic significantly. In particular, the electronic-vibrational coupling results in an 
earlier onset of the current, a quenching of the NDR effect as well as step-like substructures 
which can be associated with vibrational excitation in the B state. The earlier onset of the 
current is due to the fact that upon interaction with vibrations the resonance peaks are 
shifted to higher energies and thus enter the voltage window at lower voltage. If the voltage 
window is large enough to comprise all vibronic features, the difference between purely 
electronic and vibronic results for the current vanishes. The vibrational substructures are a 
result of the splitting of each electronic resonance into several vibronic resonances (cf. Fig. [3] 
(a)), which contribute with different weights (determined by the respective Franck-Condon 
factors) to the transmission. In contrast to the purely electronic case, the current thus 



increases in several steps. This effect is well known from previous model studies 35 
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hole injection energy (eV) 

FIG. 5: Transmission of BDET at voltages 3.94 V (solid line) and 4.32 V (dashed line) based 
on a purely electronic calculation. The two vertical lines indicate the respective lower integration 
boundaries. 

Smaller structures in the current- volt age characteristics are more clearly revealed in the 
differential conductance, dl/dV. Fig. H] shows the differential conductance of BDET ob- 
tained by numerical differentiation of the current in Fig. [3] (b). Overall, the differential 
conductance reflects the transmittance of the junction. However, in contrast to the trans- 
mission probability at zero voltage (depicted in Fig. [3] (a)), it also incorporates its voltage 
dependence and the effects of the Fermi distribution in the leads. In the purely electronic 
conductance (dashed line), this difference manifests itself in the splitting of the peak due 
to state B as well as the negative differential resistance between 4 and 4.3 Volt. At 4.5 V 
a narrow resonance due to state A can be seen. The differential conductance based on the 
vibronic calculation exhibits pronounced structures due to vibronic coupling. These include 
the shifted onset (corresponding to the 0° transitions) of state A at 3.58 V and state B at 
4.27 V as well as resonance peaks at 3.98 V corresponding the {yd)}) vibrational excitation in 
state B and at 4.32 V reflecting the (v a )o excitation in state A. The differential conductance 
based on the vibronic calculation shows, furthermore, small regions of negative differen- 
tial resistance for larger voltages, which do not exist in the purely electronic differential 
conductance. 

Finally, we consider the mechanisms responsible for the NDR effect observed in the 
results discussed above. As was mentioned above, the NDR effect is caused by the voltage 
dependence of the self energies. These result in a voltage dependence of the transmission 
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probability which in turn determine the current via the integral in Eq. (12.291) . To elucidate 
the mechanism causing NDR in the purely electronic calculation (dashed lines in Figs. [3] (b), 
H|), Fig. depicts the transmission probability for voltages of 3.94 V (corresponding to the 
local maximum in the I-V curve) and 4.32 V (the local minimum in the I-V curve). The 
peak caused by orbital B in the transmission for 3.94 V (dashed line) is much thinner than 
the corresponding peak in the curve for 4.32 V (solid line). Although the current in the 
former case corresponds to the integral with a lower boundary (see vertical line in Fig. [5]), 
the area under the '4.32 V'-transmission-curve is much smaller thus explaining the NDR in 
the current-voltage characteristic and the conductance. It should be emphasized that this 
NDR effect can only be described, if the energy and voltage dependence of the self energies 
is taken into account and will be missed within the often used wide-band approximation. 

B. Benzenediethanthiolate with a cuboid-shaped gold-cluster geometry 

To investigate the influence of the gold-cluster geometry on the transport properties of 
the molecular junction, we have studied a system where the molecule BDET is bound to 
gold-clusters with different geometry of cuboid shape (cf. Fig. [1] middle panel). To build 
this geometry we started from the optimized two-layer system described above. Next, the 
number of gold atoms in the second layer was increased from 5 to 12 and two more gold- 
(111) layers consisting of 12 gold atoms each were added, thus resulting in a metal cluster of 
38 gold atoms on each side. For better comparison with the results obtained in Sec. IIII Al 
the same value was used for the Fermi energy. The four normal modes and the molecular 
electronic states included in the calculation are very similar to those in Sec. IIII Al 

Fig. M (a) shows the zero- voltage transmission probability for this system. The comparison 
with the corresponding result for the pyramidal gold-cluster geometry (Fig. [3]) reveals overall 
a significant effect of the cluster geometry. Considering specifically the transmission based 
on the purely electronic calculation, the resonance due to state B is split into three peaks 
at -1.71 and -1.79 and -1.95 eV. This is caused by the electrode-mediated interaction of 
state B with both lower and higher lying states, which is more pronounced in the cuboid 
gold-cluster geometry. In addition, there is a small broad peak at 0.39 eV due to two sulfur 
p-orbitals, which is missing in the pyramidal geometry. The resonance due to state A in 
the zero-voltage transmission, on the other hand, is almost unaffected by the geometry of 
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the gold cluster. This is to be expected because state A is exclusively localized on the 
benzene ring. The coupling of states A and B to the contacts is of comparable size as for the 
pyramidal gold clusters: 2.17) = 4.0- 10~ 4 eV (corresponding to a lifetime of r = 1646 

fs), r BB (-1.98) = 6.5 • l(T 2 eV (r = 10 fs), and r BB (-1.71) = 1.1 ■ lO^eV (r = 6 fs). 

Considering the results of the vibronic calculations (full line in Figs. [21 E]), the most sig- 
nificant difference between the two cluster geometries is the appearance of the Op-transition 
and vibrational progressions due to state A. These were not visible in the system with the 
pyramidal gold clusters due to accidental coincidences with transmission features of other 
states. In the system with cuboid cluster geometry, the 0[]-line due to state A appears at 
-2.03 eV and the (js)q transition peaks corresponding to all four vibrational modes can be 
seen at -2.10 (a), -2.17 (b), -2.18 (c), and -2.24 (d) eV. 

Fig. [6] (b) shows the current-voltage characteristic for the BDET system with cuboid 
gold clusters. The current- volt age characteristic based on the purely electronic calculation 
(dashed line) has a similar appearance as in the system with the pyramidal cluster (cf. 
Fig. El (b)). Noticeable differences are the earlier onset and the overall larger value of the 
current as well as the appearance of two (instead of one) maxima with corresponding regions 
of NDR. The larger current is a manifestation of the overall larger molecule-lead coupling 
provided by the binding to the cuboid shaped gold cluster. The earlier onset of the current 
is caused by orbitals localized on the sulfur atoms. The appearance of two maxima is due 
to the splitting of the resonance peak of state B. We next consider the current based on the 
vibronic calculation (full line in Fig. [6] (b)). In contrast to the system with the pyramidal 
gold custer geometry, vibrational substructures are less pronounced and the current obtained 
with the vibronic calculation is in some regions significantly larger than that based on the 
electronic current. Furthermore, NDR effects are only partly quenched by the interaction 
with the vibrational degrees of freedom. 

To conclude this section, the comparison of the data for BDET bound to two different 
gold clusters shows that the geometry of the gold contact may have a significant influence on 



63, which show that the 



the transport properties. This is in accordance with results of Ref. 
transport characteristics of benzenedithiol depends on the surface structure and the lateral 
size of the electrodes. Our results show that the influence of the gold contact geometry is 
particularly pronounced for transport through molecular orbitals which have a significant 
contribution at the boundary to the contacts such as state B in the system studied here. 
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FIG. 6: (a) Total transmission probability through a BDET molecular junction with a cuboid- 
shaped gold cluster geometry at zero voltage as a function of the initial energy of the hole (relative 
to the Fermi energy). The two orbitals, denoted A and B, dominate the transmittance at the 
indicated peaks. Only the energy range with non-zero transmittance in the interval [-2.5:+2.5], 
corresponding to a voltage window of 5 V, is shown, (b) Current-voltage characteristic of BDET 
bound to cuboid-shaped gold contacts. Shown are results of calculations with (solid line) and 
without (dashed line) molecular vibrations. Only the non-zero part of the positive voltage range 
is shown. The different lines in the two panels show results of calculations with (dashed line) and 
without (solid line) vibronic coupling. 

In experiments, where it is often not straightforward to control the molecule-lead binding 
geometry, this difference may result in fluctuations of the current for different realizations 
of molecular junctions. 
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FIG. 7: Normal modes of BDBT included in the calculation. 



C. Benzenedi(butanethiolate) 



The systems studied in the previous sections are characterized by a moderate molecule- 
lead coupling, corresponding to a relatively short lifetime of the electron on the molecular 
bridge. Accordingly, vibrational effects are noticeable but not very pronounced. The cou- 
pling between the 7r-system of the molecular bridge and the gold electrodes can be reduced 
(and thus the lifetime of the electron on the bridge extended) by increasing the length of 
the alkyl-spacer group. As an example of such a system, we consider in this section elec- 
tron transport through p-benzene-di(butanethiolate) (BDBT) (cf. Fig. [TJ bottom panel). 
Thereby, the contacts were modelled with a gold cluster of pyramidal geometry. The trans- 
port properties of this system can be well described including six electronic states localized 
at the molecular bridge in the calculation. These comprise the states A and B illustrated 
in Fig. [HI which are similar as in BDET, and four sulfur p-orbitals. As in the previous 
system, from all 54 totally symmetric normal modes the four with the strongest vibronic 
coupling were included explicitly in the calculation. The modes are depicted in Fig. [71 Their 
frequencies and coupling constants are given in Table [2j 

Fig. [8] (a) shows the zero- voltage transmission probability of BDBT. In contrast to BDET, 
the transmission probability exhibits well-separated narrow resonance structures. This is 
due to the significantly smaller electronic coupling. The results based on the purely elec- 
tronic calculation exhibit two narrow resonances at -1.78 and -1.39 eV. These are due to 
states A and B, respectively, which couple to the gold electrodes with coupling strengths 
of r AA (-1.78) = 1.0 • l(T 3 eV (r = 659 fs) and r BB (-1.39) = 5.4 • l(T 4 eV (r = 1220 fs), 
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TABLE 2: Parameters of the four most important vibrational modes of BDBT between pyramidal 
gold contacts including frequencies, periods, and vibronic coupling in the electronic states A and 
B. 





uj (cm 1 ) 


T(fs) 


K (A) ( meV ) 


k (b ^> (meV) 


(a) 


566.12 


59 


69 


29 


(b) 


1198.40 


28 


52 


73 


(c) 


1232.90 


27 


114 


56 


(d) 


1676.10 


20 


132 


170 



respectively. Due to the small coupling, these states exhibit negligible level shift. Includ- 
ing vibronic coupling both peaks are shifted to the respective 0q transition and pronounced 
vibrational progressions appear. The individual peak heights correspond to the respective 
Franck-Condon factors. Table [3] gives an assignment of the most prominent resonances. 



orbital A 






mode (a) 
mode (b) 
mode (c) 
mode (d) 


-1.72 

-1.80 
-1.85 


-1.79 

-1.93 
-2.07 


orbital B 


iy)l 


MS 


mode (a) 
mode (b) 
mode (c) 
mode (d) 


-1.35 
-1.42 
-1.43 
-1.48 


-1.56 
-1.71 



TABLE 3: Assignments of the most pronounced vibronic transitions in the transmission function of 
BDBT depicted in Fig. [8] (a). The location of the features are given in eV. Blank entries correspond 
to transitions that could not be identified because of too small vibronic coupling. 

The current- volt age characteristic of BDBT is depicted in [8] (b). The result based on the 
purely electronic calculation shows two prominent steps corresponding to states A and B at 
voltages 2.76 V and 3.56 V, respectively. In addition, there is a small current for low voltages 



25 



hole injection energy (eV) 
-2.5 -2 -1.5 -1 




0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 
voltage (V) 



FIG. 8: (a) Total transmission probability through a BDBT molecular junction with a pyramidal 
gold cluster geometry at zero voltage as a function of the initial energy of the hole (relative to the 
Fermi energy). The two orbitals, denoted A and B, dominate the transmittance at the indicated 
peaks. Only the energy range with non-zero transmittance in the interval [-2.5:+2.5], corresponding 
to a voltage window of 5 V, is shown, (b) Current-voltage characteristic of BDBT bound to 
pyramidal gold contacts. Shown are results of calculations with (solid line) and without (dashed 
line) molecular vibrations. Only the non-zero part of the positive voltage range is shown. The 
different lines in the two panels show results of calculations with (dashed line) and without (solid 
line) vibronic coupling. 

and another step-like structure at 2 V. Both features are due to sulfur p-states, which cause 
a small but very broad peak in the transmission probability around the Fermi level and a 
small, narrow peak at -1 eV, which can be seen in the transmission probability at 2 V only 
(data not sown), but not in the zero- voltage transmission. As in the other model systems, 
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NDR effects caused by the same mechanism as for BDET appear at larger voltages. Upon 
including vibronic coupling, the onset of the current related to states A and B is shifted 
to lower voltage due to nuclear relaxation. Furthermore, a number of additional steplike 
structures appear which can be assigned to excitation of the different vibrational modes in 
states A and B. The most pronounced features are analyzed in Table HI The NDR following 
the last step still exists but is significantly smaller compared to the purely electronic case. 





orbital A 


orbital B 


<ya)h 
(yt)h 

(ya)h 


3.27 
3.42 
3.57 
3.57 
3.73 


2.55 
2.68 
2.85 
2.85 
2.95 
3.42 



TABLE 4: Assignments of the most pronounced vibronic transitions in the I-V curve of BDBT 
depicted in Fig. [8] (b). The location of the features are given in V. Blank entries correspond to 
transitions that could not be identified, because of too small vibronic coupling. 



IV. CONCLUSIONS 

In this paper we have studied the effect of vibrational motion on resonant charge transport 
in single molecule junctions. The methodology used to describe vibrationally coupled charge 
transport is based on a combination of first-principles electronic structure calculations to 
characterize the system and inelastic scattering theory. To apply this methodology to sys- 
tems where the transport is dominated by the occupied molecular orbitals of the junction, 
we have extended the approach to allow for hole transport. 

We have applied the methodology to molecular junctions where a benzene molecule is 
connected via alkanethiol bridges to two gold electrodes. This class of molecular junctions 
was chosen, because the change of the length of the alkyl chain allows a systematic variation 
of the coupling between the n system of the phenyl ring and the gold electrodes and thus 
a variation of the lifetime of the electron on the molecular bridge. The ratio between the 
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vibronic coupling and the electronic molecule-lead coupling determines the importance of vi- 
bronic effects in molecular junctions. The results of this study demonstrate this trend. Our 
previous work 42j had indicated that in benezendithiolate, which is characterized by rather 
strong molecule-lead coupling (and thus a very short lifetime of the electron on the molecular 
bridge), vibrational effects in resonant transport are almost negligible. Benezendibutanthio- 
late, the system with the smallest molecule-lead coupling investigated, exhibits pronounced 
vibronic effects. In particular, electronic-vibrational coupling result in a splitting of elec- 
tronic resonances into vibronic subresonances in the transmission probability. In the current- 
voltage characteristic, vibronic effects manifest themselves in steplike structures, which can 
be associated with vibrational states in the molecular cation. Benezendiethanethiolate is 
an intermediate case, which shows noticeable though not very pronounced vibronic effects. 
The study also shows that electronic-vibrational coupling may result in a quenching of neg- 
ative differential resistance effects, which are caused by the voltage dependence of the self 
energies, and a significantly altered overall magnitude of the current. 
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